#read in data
#must have the following file path in working directory: data/kernfields/kernYEAR/kernYEAR.shp

#create list of years
kern_field_layers <- c(2016:2019)

#loop through read in process
for (i in kern_field_layers) {
  filename <- paste0("kernfields_", i)
  wd <- paste0("data/kernfields/kern", i, "/kern", i, ".shp")
  assign(filename, read_sf(wd))
}

#read in kern county shapefile
kern_county <- read_sf(here("data", "Kern_County_Boundary", "Kern_County_Boundary.shp")) %>% 
  st_transform(crs = crs(kernfields_2016))
#clean names, make sure each polygon has same number of rows

kern_16_clean <- kernfields_2016 %>% 
  clean_names() %>% 
  mutate(symbol = NA) %>% 
  mutate(comm_code = NA) %>% 
  mutate(year = 2016)

kern_17_clean <- kernfields_2017 %>% 
  clean_names() %>% 
  mutate(year = 2017)

kern_18_clean <- kernfields_2018 %>% 
  clean_names() %>% 
  mutate(year = 2018)

kern_19_clean <- kernfields_2019 %>% 
  clean_names()  %>% 
  mutate(year = 2019)


#bind all rows,dissolve (union) boundaries
kern_binded <- rbind(kern_16_clean, kern_17_clean, kern_18_clean, kern_19_clean) %>% 
  st_union()

#plot
plot(st_geometry(kern_binded))

#write to sf
#st_write(kern_binded, here("data", "shapefiles_written", "unioned_fields_20162019"), driver = "ESRI Shapefile")
#crop

#select region (based off of farmer)
kern_subset <- kern_19_clean %>% 
  filter(agent == "McMANUS/WILSON,MICHELE/AARIN")

plot(kern_subset["p_status"])

#crop
kern_19_crop <- st_crop(kern_19_clean, kern_subset)
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries
plot(kern_19_crop["p_status"])

kern_18_crop <- st_crop(kern_18_clean, kern_subset)
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries
kern_county_crop <- st_crop(kern_county, kern_subset)
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries
plot(kern_county_crop["FID"])

# Convert to tibble for a dplyr full join on the geometry column.
kern1819_intersect <- full_join(as_tibble(kern_19_crop), as_tibble(kern_18_crop), by = "geometry")

kern_1819_sf <- st_as_sf(kern1819_intersect)

tmap_mode("view")
## tmap mode set to interactive viewing
tm_shape(kern_1819_sf) +
  tm_polygons(c("pmt_site.x", "pmt_site.y"), legend.show = FALSE)
## Warning: Number of levels of the variable "pmt_site.x" is 459, which is
## larger than max.categories (which is 30), so levels are combined. Set
## tmap_options(max.categories = 459) in the layer function to show all levels.
## Warning: Number of levels of the variable "pmt_site.y" is 478, which is
## larger than max.categories (which is 30), so levels are combined. Set
## tmap_options(max.categories = 478) in the layer function to show all levels.
tm_shape(kern_18_crop) +
  tm_polygons("pmt_site", legend.show = FALSE)
## Warning: Number of levels of the variable "pmt_site" is 478, which is
## larger than max.categories (which is 30), so levels are combined. Set
## tmap_options(max.categories = 478) in the layer function to show all levels.
tm_shape(kern_19_crop) +
  tm_polygons("pmt_site", legend.show = FALSE)
## Warning: Number of levels of the variable "pmt_site" is 459, which is
## larger than max.categories (which is 30), so levels are combined. Set
## tmap_options(max.categories = 459) in the layer function to show all levels.
plot(st_geometry(kern_19_crop), main = "2019")

plot(st_geometry(kern_18_crop), main = "2018")

plot(st_geometry(kern_1819_sf), main = "Intersect")

library(stampr)
library(sp)
library(smoothr)
## 
## Attaching package: 'smoothr'
## The following object is masked from 'package:stats':
## 
##     smooth
#rbind inputs, create ID column
croped_fields_binded <- rbind(kern_18_crop, kern_19_crop) %>% 
  st_make_valid() %>% 
  rownames_to_column(var="ID")

#conver to spatial
cropped_sp <- as(croped_fields_binded, "Spatial")

#use the stamp function to assess space-time change
stamp_test <- stamp.multichange(cropped_sp, changeByField = T, changeField = "year", dc = 0, distance = TRUE, direction = F)

#drop geometry "crumbs"
#stamp_test_crumbs <- drop_crumbs(stamp_test, threshold = 1, drop_empty = T)

#convert back to sf
stamp_test_sf <- st_as_sf(stamp_test) %>% 
  st_make_valid()

head(stamp_test_sf)
## Simple feature collection with 6 features and 11 fields
## Geometry type: GEOMETRY
## Dimension:     XY
## Bounding box:  xmin: 6185027 ymin: 2366668 xmax: 6187962 ymax: 2386572
## Projected CRS: NAD83 / California zone 5 (ftUS)
##        ID1 ID2 LEV1 LEV2 LEV3     LEV4 GROUP       AREA   CENDIST TGROUP
## 1-1882   1  NA DISA CONT CONT DIVISION     1   37738.53  524.4832      1
## 1-1883   1 622 STBL STBL STBL DIVISION     1 2870084.22 1033.4926      1
## 1-1884   1 623 STBL STBL STBL DIVISION     1 2291794.97 1285.7153      1
## 1-1885  NA 622 GENR EXPN EXPN DIVISION     1   34751.97 1239.5939      1
## 1-1886  NA 623 GENR EXPN EXPN DIVISION     1   22019.28 1304.4524      1
## 1-1887  13  NA DISA CONT CONT     BOTH     2       0.00  252.6785      1
##        STGROUP                       geometry
## 1-1882      11 MULTIPOLYGON (((6187962 238...
## 1-1883      11 MULTIPOLYGON (((6187962 238...
## 1-1884      11 MULTIPOLYGON (((6187962 238...
## 1-1885      11 MULTIPOLYGON (((6186808 238...
## 1-1886      11 MULTIPOLYGON (((6187926 238...
## 1-1887      12 LINESTRING (6185032 2367243...
#task is somewhat computationally intensive...write to shapefile after

#st_write(stamp_test_sf, here("data", "shapefiles_written", "stamp_test.shp"), driver = "ESRI Shapefile")
#THIS APPEARS TO WORK...NEED TO DROP CRUMBS
plot(st_geometry(kern_19_crop), main = "2019")

plot(st_geometry(kern_18_crop), main = "2018")

plot(st_geometry(stamp_test_sf), main = "Stamp_test")

library(tmap)

tmap_mode("view")
## tmap mode set to interactive viewing
tm_shape(stamp_test_sf) +
tm_polygons("GROUP")

STAMP interp

ID1 Polygon ID from T1 polygons; NA if it did not exist, ID2 Polygon ID from T2 polygons; NA if it did not exist, LEV1 Level 1 STAMP designation, LEV2 Level 2 STAMP designation, LEV3 Level 3 STAMP designation, LEV4 Level 4 STAMP designation, GROUP Group ID signifying group membership, AREA Polygon area in appropriate areal units, – (optional) Additional columns from directional analysis if direction = TRUE, – (optional) Additional columns from distance analysis if distance = TRUE,

STAMP events are reported at four levels of increasing complexity: LEV1 – disappearance (DISA), stable (STBL), and generation (GENA); LEV2 – disappearance (DISA), contraction (CONT), stable (STBL), expansion (EXPN), and generation (GENR); LEV3 – disappearance (DISA), T1 displacement (DISP1), convergence (CONV), concentration (CONC), contraction (CONT), stable (STBL), expansion (EXP), fragmentation (FRAG), divergence (DIV), T2 displacement (DISP2), and generation (GENR); LEV4 – LEV4 is different from other levels. It is used to identify those groups where union (UNION), division (DIVISION), and both union and division (BOTH) events occur.